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Abstract 

The objective of this research was to analytically 
and experimentally study the capabilities of piezo- 
electric plate actuators for suppressing flutter. Piezo- 
electric materials are characterized by their ability 
to produce voltage when subjected to a mechanical 
strain. The converse piezoelectric effect can be uti- 
lized to actuate a structure by applying a voltage. 
For this investigation , a two- degree-of- freedom wind 
tunnel model was designed, analyzed, and tested. 
The model consisted of a rigid wing and a flexible 
mount system that permitted a translational and a 
rotational degree of freedom. The model was de- 
signed such that flutter was encountered within the 
testing envelope of the wind tunnel. Actuators made 
of piezoelectric material were affixed to leaf springs 
of the mount system. Command signals, applied 
to the piezoelectric actuators, exerted control over 
the damping and stiffness properties. A mathemati- 
cal aeroservoelastic model was constructed by using 
finite element methods, laminated plate theory, and 
aeroelastic analysis tools. Plant characteristics were 
determined from this model and verified by open loop 
experimental tests. A flutter suppression control law 
was designed and implemented on a digital control 
computer. Closed loop flutter testing was conducted. 
The experimental results represent the first time that 
adaptive materials have been used to actively sup- 
press flutter. They demonstrate that small, care- 
fully placed actuating plates can be used effectively 
to control aeroelastic response. 

Introduction 

Flutter, an interaction between the structural dy- 
namics and the aerodynamics that results in diver- 
gent and destructive oscillations of motion, has been 
observed and documented on aircraft since the era of 
controlled flight began (ref. 1). Historically, passive 
solutions such as increasing structural stiffness, mass 
balancing, or modifying geometry have been used 
to prevent this hazardous phenomenon. These ap- 
proaches result in increased weight and cost, and de- 
creased performance. During the past 20 years, there 
has been considerable research to develop active flut- 
ter suppression concepts that use conventional lead- 
ing and trailing edge aerodynamic control surfaces 
(refs. 1-5). An active control approach eliminates 
most of the weight and performance penalties associ- 
ated with the passive approach and additionally pro- 
vides a flexibility so that the control law can vary 
with configuration or flight condition. Active flut- 
ter suppression is not a common practice in today’s 
commercial or military aircraft because of several 


concerns. Flutter is generally of a catastrophic na- 
ture; therefore a failure of the system could affect 
flight safety. As a result, system redundancy, reli- 
ability, and maintainability are critical issues to be 
addressed. To a lesser extent, the control surface au- 
thority available to maneuver the aircraft with the 
simultaneous implementation of active flutter sup- 
pression is also a concern. To alleviate these con- 
cerns, alternatives to utilizing the aerodynamic con- 
trol surfaces for active flutter suppression are being 
studied. 

The use of secondary controllers made of adap- 
tive material is one such concept. There arc several 
classifications of adaptive materials, including piezo- 
electrics, electrostrictors, shape memory alloys, and 
magnetostrictors. A detailed account of the prop- 
erties, benefits, and drawbacks of each type can be 
found in reference 6. The present study focused ex- 
clusively on the use of piezoelectric materials, which 
were chosen based on their wide control bandwidth, 
material properties, and availability of results from 
previous tests. Piezoelectricity is the ability of a 
material to develop an electrical charge when sub- 
jected to a mechanical strain. The converse piezo- 
electric effect, the development of mechanical strain 
when subjected to an electrical field, can be utilized 
to actuate a structure. A local strain is produced 
in the structure, which induces forces and moments. 
By judicious arrangement of piezoelectric plates, the 
correct reaction of the structure required to inhibit 
flutter can be produced. Many research efforts have 
utilized adaptive plate actuators for various applica- 
tions (refs. 6-17), and they have been recently used 
in an analytical flutter suppression study (ref. 8). Re- 
sults available from aeroelastic applications of piezo- 
ceramics are very limited. Static aeroelasticity has 
been the subject of investigations by Ehlers and 
Weisshaar (refs. 6, 15, and 16). They conducted an- 
alytical studies on laminated composite wings with 
embedded actuators, looking at pure torsional and 
bending deformations. They reported that through 
feedback to embedded adaptive material layers, the 
divergence speed is altered, implying also that lift 
effectiveness is influenced. The augmentation or re- 
placement of conventional aerodynamic control sur- 
faces with strain actuation for aeroelastic control has 
been the focus of an analytical investigation of a typ- 
ical section by Lazarus, Crawley, and Lin (ref. 17). 
They found that strain actuation by means of piezo- 
electric elements may provide a viable and effective 
alternative to articulated control surfaces for control- 
ling aeroelastic response. Investigation of flutter sup- 
pression for lifting surfaces and panels has been done 
by Scott (ref. 8). That analytical study considered 



controlling flutter at supersonic speeds by using full 
state feedback. 

The purpose of the present study was to investi- 
gate flutter suppression by using piezoelectric plates 
as actuators. Results from analyses and experiments 
demonstrating this technology will be presented. A 
discussion of piezoelectric materials and some appli- 
cations will be presented first. Following this, details 
of the experimental configuration and testing facil- 
ities are given. Development and implementation 
of the aeroservoelastic equations of motion are pre- 
sented. Results of analytical studies based on these 
open loop aeroservoelastic equations, the design of 
the control law, and results from closed loop studies 
are next presented. Experimental results and com- 
parison with analytical predictions follow. System 
identification test results are discussed as well as the 
open and closed loop flutter results. 

Symbols 


A 

area 

c 

chord 

^15)^24 

electromechanical coupling 
coefficient, shearing response to 
voltage applied perpendicular to 
the poling axis 

^31 , ^32 

electromechanical coupling 
coefficient, in- plane response to 
voltage in poled direction 

^33 

electromechanical coupling 
coefficient, response in the 
direction of applied voltage 
corresponding to the poling axis 

e 3 

applied voltage per thickness, in 
the 3-direction 

E 

Young’s modulus 

f 

natural frequency in cycles per 
second (hertz) 

G.F. 

gauge factor 

G 12 

shear modulus 

9 

structural damping 

lyv 

inertia about the pitch axis 

k 

stiffness of a single degree of 
freedom 

m 

mass 

™ 0 

mass per unit area 


n 

P 

Q 

Qxii Qyii Qzi 


T 

t 

t 

U 

u 


U , V, w 

V 

V in ax 

V 

x,y,z 

x n 


x 1 

z 


7 

7 xy 

S 

e 

£xx 

£ yy 

C 

K 

A 

V 

p 

a 


number of elapsed cycles 

transfer function 

dynamic pressure 

generalized coordinates corre- 
sponding to mode shape displace- 
ment vectors in the x-, y-, and 
2 -directions, respectively 

Laplace variable 

kinetic energy 

thickness of the piezoelectric 
plate 

time 

potential energy 

control command in the state 
space equations 

component displacements in 
the x-, y-, and 2 -directions, 
respectively 

applied voltage 

maximum voltage that may be 
applied to the actuator 

velocity 

component directions 

amplitude of the nth cycle 

location of the Oth point within 
the x, y, 2 frame 

amplitude of the first cycle 

distance from the midplane of the 
plunge spring tine to the bonded 
face of the piezoelectric plate 

piezoelectric actuation constant 
shear strain in the x-y plane 
logarithmic decrement 
strain 

normal strain in the x-direction 

normal strain in the y-direction 

damping ratio 

curvature 

eigenvalue 

Poisson’s ratio 

density 

stress 
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G x 

normal stress in the x-direction 

M 

generalized mass matrix 

(Jy 

T X y 

normal stress in the y-direction 

M 

physical mass matrix (not 

shear stress in the x-y plane 


generalized) 

UJ 

natural frequency 

M 

generalized mass matrix that 

Matrices 

A 

A 

and vectors: 

state matrix 

extensional stiffness matrix 

m a 

includes the structural and 
aerodynamic second order terms 

moment due to strain actuation 

force due to strain actuation 


A2 


B 

B 

C 

D 

D 

D 

D 

d 

E 

Pi 

F 2 

f 3 

G 

K 

K 

K 


constant terms of the rational 
function approximation to the 
generalized aerodynamic forces 

coefficients of the first derivative 
terms of the rational function 
approximation to the generalized 
aerodynamic forces 

coefficients of the second deriva- 
tive terms of the rational func- 
tion approximation to the gener- 
alized aerodynamic forces 

state space control matrix 

coupling stiffness matrix 

state space output matrix 

generalized damping matrix 

physical damping matrix (not 
generalized) 

generalized damping matrix 
that includes the structural and 
aerodynamic first order terms 

bending stiffness matrix 

electromechanical coupling 
matrix 

applied voltage per thickness 

potential energy term 

strain actuator term 

kinetic energy term 

stiffness matrix 

generalized stiffness matrix 

physical stiffness matrix (not 
generalized) 

generalized stiffness matrix 
that includes the structural and 
aerodynamic zeroth order terms 


Q 

Q 

q 

R 

S 

u 

u 

X 

Xfi,Xf 2 

y 

P 

T 

£ 

€ 

K 

A 

(7 

$UU 

^uy 

* 

^xu^yi^zi 

Subscripts: 

a 

f 


aerodynamics 

rational function approximations 
to the aerodynamic generalized 
forces 

vector of generalized coordinates 
electric flux density 
compliance matrix 
control vector 

displacement vector with compo- 
nents u , u, and w 

state vector 

states associated with the flexible 
modes and their derivatives 

output vector 

permittivity matrix 

actuating matrix 

strain vector 

impermittivity matrix 

curvature vector 

actuation strain vector 

stress vector 

autospectrum of input 

cross-spectrum of input with 
output 

mode shape displacement matrix 

mode shape displacement vectors 
for node point i in the x-, y-, 
and z-directions, respectively 


associated with the aerodynamics 

associated with the flexible 
modes 
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ij 

ith and jth node point, mode, or 
element 

L , M , and N 

limits in the x-, y- and 
^-directions, respectively 

V 

associated with the controls or 
piezoelectrics 

pi 

associated with the plunge degree 
of freedom 

s 

associated with the structure 

x,y,z 

quantities in the x-, y-, and 
2-directions, respectively 

0 

associated with the pitching 
degree of freedom 

A 

associated with the strain 
actuation 

0 

associated with the midplane 

Operators: 


V 

derivative operator matrix 

dA 

dummy variable of integration 
over the area 

dV 

dummy variable of integration 
over the volume 

{ } 

vector 

n 

matrix 

[ ) T 

transpose 

0 

derivative with respect to time 


Piezoelectric Materials 

A material which, when subjected to a mechan- 
ical load, accumulates an electric charge is said to 
have piezoelectric properties. Many naturally occur- 
ring crystals have piezoelectric capabilities, but their 
nonuniform properties encouraged research into man- 
ufacturing materials that would produce electrome- 
chanical coupling (ref. 18 ). Certain polymers and 
ceramics consist of crystalline subdomains that are 
bipolar in nature. Piezoelectric properties can be 
induced in these materials by applying a large elec- 
trical field across them. This induces an orienta- 
tion of crystalline subdomains such that the posi- 
tive and negative poles of the individual domains are 
aligned with the applied voltage field, denoted the 
3 -direction. Materials that require this poling belong 
to a subcategory termed ferroelectrics. 


The poling process of ferroelectric materials re- 
orients the dipoles such that there is a net polariza- 
tion along the axis of applied voltage. The orien- 
tation remains after the inducing field is removed. 
Using these materials requires a voltage to be subse- 
quently applied through electrodes on opposite faces 
of the material. The most common configuration 
(fig. 1) is to place the electrodes on the faces par- 
allel to the poled axis and to apply the voltage in the 
same direction as the original inducing field. The ma- 
terial deforms both through the thickness, denoted 
the ds 3 effect (fig. 2), and in the in-plane directions, 
denoted the c/31 or c/32 effect (fig. 3 ). Applying a 
field oriented in one direction induces in-plane ex- 
pansion; applying it in the other direction induces 
in-plane contraction. To define the electromechani- 
cal effects, the first subscript denotes the direction 
of the applied voltage and the second subscript de- 
notes the direction of the deformation. For trans- 
versely isotropic materials, there is no distinction be- 
tween vectors lying in any plane perpendicular to the 
poling axis. To complete the description of the de- 
formations achievable with these materials, figure 4 
shows the electrodes placed on faces parallel to the 
poling axis (i.e., in the 1 -direction). This induces a 
shearing strain within the piezoelectric, as the posi- 
tively poled side of the piezoceramic strains toward 
the negatively charged electrode and the negatively 
poled side strains oppositely. 



Figure 1. Electrode placement on piezoelectric plate. 


Choosing the proper piezoelectric material to use 
for a given application is based on stiffness prop- 
erties, flexibility, electromechanical coupling coeffi- 
cients, and limits on applied voltage. The ability of a 
piezoelectric material to actuate a structure is a func- 
tion of its stiffness, the limit on the voltage that can 
be applied across it, and the electromechanical cou- 
pling coefficients. Polymers have high voltage limits, 
yet they have low stiffness and low electromechanical 
coupling coefficients. Ceramics, on the other hand, 
are much stiffer and have large coupling coefficients 
and are thus better suited for actuator applications. 
Lead zirconate titanate (PZT), a piezoceramic, was 
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Figure 2. Thickening effect, (f/.w effect). 



Figure 3. Lengthening effect (f/;u , dyi effect). 



Figure 4. Shearing effect {d\^, di\ effect). 


chosen for this investigation. The material properties 
(Product Information Catalog, PiezoSystems Solid 
State Motion Technologies) of the most concern in 
this study are listed in table I. 


Table I. Material Properties for Lead 
Zirconate Titanate (PZT) 


Property 

Value 

d:u 

dx\ 

Maximum voltage 
Young’s modulus 
Poisson's ratio 
Density 

166 x 10“ 12 m/V 
-360 x Kr 12 m/V 
0.7 x 10 K V/m 
63 x 10 9 N/m 2 
0.33 

7.65 x 10'* kg/m' ! 


Piezoelectric plates can be configured in different 
ways to accentuate the displacements or forces being 
generated. The in-plane expansion and contraction 
of adaptive materials may be utilized by bonding 
actuating plates to either side of a center shim (fig. 5). 
One plate is expanded and one is contracted; the 
net result is a bending displacement much greater 
than the length deformation of either of the two 
layers. This configuration, which takes advantage of 
the Poisson-like d%\ effect, is referred to as a bimorph 


Undeformed 



Figure 5. Bimorph bender configuration of piezoelectric plate 

elements. For the voltage polarity shown, the top layer 

expands in the in-plane direction and the bottom layer 

contracts in the in- plane direction. 

bender element. It will serve as the primary actuator 
mechanism for the investigation described herein. 

The material properties of the piezoelectric ma- 
terials have been treated in this work with lin- 
ear relationships. This assumption is valid for low 
applied voltages and small deformations. Nonlinear- 
ities and nonidealities of these materials have been 
well-documented in references G, 7, and 9. Several 
nonideal properties that have been found to have sig- 
nificance arc the amplitude dependence of the field- 
strain relationship, creep, variations with mechani- 
cal strain, and depoling. These issues will not be 
addressed in detail here; however, efforts have been 
exerted during this investigation to avoid these non- 
linear regions. The amplitude of the control signal 
voltage was low, thus depoling was avoided and the 
linear field-strain relation was maintained. The fre- 
quency of concern was approximately 10 Hz, thus 
creep, a low-frequency phenomenon, was avoided. 
The plates were also placed on a region where the me- 
chanical strain at rest was very small, avoiding non- 
linearities associated with applied mechanical strain. 

Piezoelectric materials can be used in sensing and 
actuating applications. Piezoelectric devices used as 
sensors emit voltages when subjected to a mechani- 
cal load. Because this study utilized the piezoelectric 
plates only as actuators, sensor applications will not 
be discussed further. In an actuating application, 
the converse piezoelectric effect is utilized as the ac- 
tuators deform in response to a control signal or ap- 
plied voltage. The mention of actuators suggests the 
use of hydraulics, pistons, motors, etc. An actuator. 
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however, is a device that moves or causes action of 
something else. The use of adaptive materials in this 
manner has lassoed engineering interest from various 
areas. The following section provides an overview of 
several investigations that have been conducted. 

In the area of rotorcraft, two distinctly differ- 
ent actuator configurations have been examined for 
higher harmonic control (refs. 9 and 10). The first 
used directionally attached plates to torsionally acti- 
vate blade sections and actuate a trailing edge flap. 
The magnitude of flapping vibrations was signifi- 
cantly reduced with active controls. The second uti- 
lized a push-pull configuration of bender elements. 
Another actuating application, also detailed in refer- 
ence 10, is the active damping of truss members for 
large space structure applications. This study used 
commercially available actuators that utilize the c/33 
effect (the expansion direction coincides with the di- 
rection of polarization) to limit the vibration ampli- 
tude and settling time of transients induced by dy- 
namic perturbations to the structure, such as crew 
motion. In the acoustics field, recent work (ref. 12) 
has focused on reducing cabin noise through de- 
structive interference produced by distributed piezo- 
electric plates. Separate finite impulse response 
filters were constructed to control an acoustic res- 
onance and a structural mode occurring 25 percent 
above the acoustic resonance. Reference 13 describes 
the use of piezoelectric plates in a biinorph configu- 
ration on an aluminum beam in conjunction with an 
adaptive controller to attenuate vibrations with fre- 
quencies above 300 Hz. Reference 14 details exper- 
iments and analyses of a composite beam with dis- 
tributed embedded actuators controlling structural 
modes from 11 to 150 Hz. Through feedback of ve- 
locity, structural damping was increased by an order 
of magnitude. 

Experimental Apparatus 

The hardware involved in this wind tunnel test 
is described in four sections: the wind tunnel, the 
test article, the digital controller, and additional 
instrumentation. 

Wind Tunnel 

The piezoelectric flutter suppression model was 
tested in the Flutter Research and Experiment De- 
vice (FRED) at Langley Research Center, shown in 
figure 6. The FRED is an open circuit tabletop 
wind tunnel with a maximum operating velocity of 
85 mi/hr (approximately 1500 in/sec). The test sec- 
tion is 6 in. by 6 in. and is constructed of acrylic 
glass for model viewing. The air is pulled through 


the tunnel by a 2-hp motor and smoothed by a single 
honeycomb screen at the beginning of the contraction 
duct. Models are mounted from the removable ceil- 
ing of the test section. Figure 7 shows the test article 
suspended from the test section ceiling. 

Wind 



Test section: 6 in. by 6 in. 

Models suspended from ceiling 
Limits: Velocity-1500 in/sec (85 mi/hr) 

Dynamic pressure-0.13 Ib/irr 

Figure 6. Schematic of Flutter Research and Experiment 

Device (FRED) wind tunnel. 

Test Article 

The piezoelectric flutter suppression wind tunnel 
model (fig. 7) consists of four integral components: 
a flexible mount system, a rigid wing, piezoelectric 
plate actuators, and a strain gauge bridge. The test 
article was designed through an iterative procedure 
involving parametric studies, the details of which are 
presented in appendix A. There were three driving 
factors in the design: the model had to flutter within 
the wind tunnel envelope, had to fit inside the wind 
tunnel with certain margins of safety, and had to 
have flat surfaces on which piezoelectric actuating 
plates could be mounted. It was decided a priori 
that a flexible mount system would reside outside the 
wind tunnel and provide the degrees of freedom for 
a rigid wing. The physical system has two distinct 
structural modes. When no mass coupling is present, 
the lower frequency mode is a plunge mode (out-of- 
plane translation of the wing). The higher frequency 
mode is a rotation about the wing pitch axis. The 
wing was designed such that the modes would be 
coupled by virtue of the wing mass distribution. 

The mount system, shown in figures 8 and 9, 
suspended the wing by two pins through slots in the 
test section ceiling and provided the plunging and 
pitching freedoms by virtue of separate spring tine 
mechanisms. The plunge mechanism consists of two 
spring steel plates or tines separated by 0.75 in. and 
clamped at both ends to maintain this distance. This 
provides the pure plunging motion of a beam with 
guided boundary conditions instead of the flapping 
motion associated with a cantilevered beam. The 
pitch mechanism is a single spring tine connected to 
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the wing at the leading edge and at the 0.2353-chord 
location, where there is a bearing-like mechanism 
that allows for free rotation. This configuration 
provides the wing with pitch stiffness and a pitch 
axis. The two mechanisms are joined together as 
shown in figure 10. The forward end of the pitch 
mechanism is fixed relative to the plunge springs by 
mounting the pitch pivot pin to the lower clamping 
block of the plunge mechanism. 


Plunge 



balsa wood extension, and an aluminum mass bal- 
last. The primary wing structure is formed from 
Vs-in-thick isotropic aluminum with a diamond cross 
section and blunted leading and trailing edges and 
midchord section. It has a chord of 2 in., with the 
pitch pivot at the midchord. The balsa wood exten- 
sion overlays the aft half of the primary structure 
and extends the chord length to 4.25 in. The trail- 
ing edge of this section was coated with aluminum to 
provide a mass ballast. The mass of the entire wing 
is 0.090 lbm and has an inertia about the pitch axis 
of 0.134 lbm-in 2 . All three sections extend the full 
span of the wing, which is 4 in. Table II gives the 
measured mass and location of the center of grav- 
ity for each portion of the wing. Based on measured 
dimensions, mass, and distance to the pitch point, 
inertias for the component parts were calculated; the 
results are given in table III, where 

leg - ^ (Mass) (Width) 2 (1) 

/pivot = (Mass) (Distance) 2 + I cg (2) 

Table II. Measured Mass and Center-of-Gravity Locations 
for Wing Components 


Figure 9. Mount system close-up: interconnection mechanism 
for the plunge and pitch degrees of freedom. 


Adjusts up and 
down to change 
effective length of | 
plunge spring 


Piezoelectric 

actuator 



Velocity 


Figure 10. Wind tunnel mount system (not to scale). 

The wing, depicted in figures 11 and 12, consists 
of three sections: an aluminum primary structure, a 



Mass, 

lbm 

Center- of- gr a vi t y 
distance aft of 
pivot point, in. 

Primary wing structure 

0.072 

0.0 

(includes pivot 



mechanism) 



Balsa extension 

.011 

1.875 

(includes adhesives) 



Aluminum mass ballast 

.007 

3.0 

Total 

0.090 



Table III. Inertia Calculations for Wing Components 



M ass . 
lbm 

Width, 

in. 

Distance, 

in. 

/eg, 

lbm-in 2 

f pivot. ? 

lbm-in 2 

Primary 

wing 

0.072 

2.0 

0.0 

0.022 

0.022 

Wing 

extension 

.011 

3.25 

1.875 

.0097 

.0487 

Mass 

ballast 

.007 

.5 

3.0 

.00014 

.06314 
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Figure 1 1 . Photograph of wing. 
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Balsa wood 
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edge 


Aluminum 

mass 

ballast 


(a) Cross-sectional view. 


(b) Planforrn view. 


Figure 12. Schematic of wing. 
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Piezoelectric ceramic plates were installed near 
the root of one of the plunge spring tines to actuate 
the test article. Two plates, 1.5 in. long and 1 in. 
wide, were bonded to opposing sides of the plunge 
spring tine, with their positive poles both oriented 
toward the steel, to form a bimorph-type actuator. 
The 0.0075-in-thick plates were electrically isolated 
from the steel by the bonding layers. Small copper 
tabs affixed beneath the plates during the bonding 
process served as the means of applying voltages to 
the bonded-side electrodes (fig. 13). 

Steel plunge 



-V 


Figure 13. Piezoelectric actuator attachment to plunge spring 

tine. 

A strain gauge bridge was mounted near the base 
of one of the spring tines, with two gauges on either 
side of the tine. The gauges can be seen in the 
photograph of figure 7. The gauges, which had 
an overall gauge factor of 2.075, were configured to 
measure the cantilever bending strain. The bridge 
was powered by a ±5-V power supply. The strain was 
computed by taking the ratio of the output voltage 
to the input voltage and dividing by the gauge factor. 


The output voltage from the strain gauge was ampli- 
fied by 100 before being sent to the digital computer. 

Digital Controller 

The control law is implemented by using a 
personal computer, with an 80386 processor and 
an 80387 coprocessor running a real-time Unix op- 
erating system. The control laws are programmed 
in the C-language and use floating point arithmetic 
for all control law calculations. The data acquisi- 
tion system uses 12-bit analog-to-digital converters 
(ref. 19). 

Additional Instrumentation 

The wind tunnel and model had three sensor sys- 
tems: a strain gauge bridge, a linear accelerometer, 
and a hot wire anemometer. The strain gauge bridge 


utilized to generate the feedback measurement signal 
has been discussed previously. 

An Endevco piezoelectric accelerometer was used 
in system identification experiments. It was powered 
by an external 4-mA current source. The output was 
calibrated at 9.98 m V/g. The accelerometer served 
as a roving sensor, being placed where applicable for 
different experiments. During zero-airspeed testing, 
it was located on the wing; however, during flutter 
testing, it was installed on the clamping block. 

Both the strain and the acceleration were ampli- 
fied by 100 before being sent to the analog-to-digital 
converters. The output voltage of the controller was 
sent to an operational amplifier having a gain of 25 
and a limit on the output voltage of 80 V. This lim- 
ited the usable range of output values from the con- 
troller to ±3.2 V. The amplified voltage was then 
applied across each of the piezoelectric elements. An 
Apex Microtechnology P83A operational amplifier 
was used to boost the input voltage to the piezo- 
electric actuators. The signal source or input voltage 
applied to the piezoelectric plates was amplified by 
a factor of 25, with a limit on the output voltage 
equal to the power supply voltage, which in this ex- 
periment was ±80 V. A block diagram of the closed 
loop system with active feedback is presented in fig- 
ure 14. Only the signals employed in the feedback 
scheme are shown. 

A Kurz 443M air velocity meter provided visual 
readouts of the test section airspeed. This is a hot 
film anemometer with an analog display in meters 
per second. The probe was inserted into the flow 
just behind the model in the test section. Thus, in 
order to accurately measure the velocity, the model 
must be moved to the stops to eliminate blockage and 
the influence of wing oscillations on the reading. 

Analytical Modeling 

The analytical model was developed incorporat- 
ing the aeroservoelastic equations of motion, a model 
of the control computer dynamics, and experimen- 
tally determined correction factors. The following 
discussion of the analytical modeling first addresses 
the theoretical development of the equations of 
motion for an elect romechanically actuated aero- 
elastic structure. Following this, the procedures used 
in generating and assembling each of the pieces of 
the open loop plant are presented. Details of com- 
puter software utilized for implementing the analyti- 
cal model are discussed. Modeling of the control com- 
puter dynamics follows, while the final section details 
scaling the model for units and amplifiers along 
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Figure 14. Block diagram of closed loop system with active feedback 


with incorporating the experimentally determined 
correction factors. 

Theoretical Development of the 

Aeroservoelastic Equations of Motion 

The aeroelastic equations of motion were de- 
veloped with classical techniques available to the 
aeroelastic community; the control forces from the 
piezoelectric elements were developed with laminated 
plate theory in conjunction with Lagrange’s energy 
method. 

Aeroservoelastic equations of motion based on La- 
grange’s energy equations represent a summation of 
forces that include the inertial, dissipation, and in- 
ternal restoring forces, the reduced- frequency depen- 
dent aerodynamic forces due to the structural mo- 
tions, and the control forces. Generalized coordi- 
nates were defined as the modal contributions to the 
displacements. An in-depth development of the gen- 
eral aeroservoelastic equations of motion is presented 
in appendix B. To highlight the difference between 
an electromechanically actuated system and an aero- 
dynamically actuated system, a brief discussion of 
the development follows. 

Solution of Lagrange’s equations requires expres- 
sions for kinetic and potential energy. The kinetic en- 
ergy is defined by a volume integral, which when in- 
tegrated over a uniform thickness laminate becomes 
an area integral, premultiplied and postmultiplied by 
the derivative of the generalized coordinates. The 


area integral is the structural mass matrix. The 
potential energy is defined by a volume integral of 
the strain and stress vectors multiplied. The con- 
stitutive relationship that defines the dependence 
of strain on stress for piezoelectric materials is de- 
rived in appendix C. The modeling of piezoelectric 
systems requires consideration of both mechanical 
and electrical behavior. Coupling between mechani- 
cal stresses and electrical fields is analytically repre- 
sented by constitutive relationships that contain both 
the electrical quantities and the mechanical quanti- 
ties. This relationship is utilized in expressing the 
potential energy; a portion of the potential energy 
expression is the structural stiffness. The presence of 
the strain due to electrical energy generates a term 
that is not traditionally observed in the mechanical 
potential energy expression. This is the piezoelectric 
actuation matrix, the development of which is pre- 
sented in appendix D. The calculation of this matrix 
is easily performed utilizing the discretized structural 
model. Derivatives of the expressions for kinetic and 
potential energy are taken as appropriate and imple- 
mented in the Lagrange equations of motion. The 
state vector is defined as the generalized coordinate 
vector and its derivative; the second order equations 
are thus transformed into first order state space form. 

Aeroservoelastic Model Construction 

A variety of techniques and software packages 
were employed to develop the various elements of 
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the equations of motion. The structural charac- 
teristics were calculated by utilizing a finite ele- 
ment method and normal modes analysis for de- 
termining the natural frequencies and mode shapes. 
Generalized aerodynamic forces were calculated from 
doublet-lattice lifting surface theory applied to the 
discretized wing. To implement the aerodynamic 
forces in state space equations, rational function ap- 
proximations were made to the tabular output of the 
doublet-lattice code. The control distribution ma- 
trix B and the output or measurement matrix C were 
computed with finite difference techniques. 

Finite element model . A finite element model 
of the wind tunnel wing and the two-degree-of- 
freedom mount system (fig. 15) was constructed and 
analyzed with MSC NASTRAN (ref. 20). The model 
developed represented the primary wing with solid 
elements; the wing extension and mass ballast were 
represented with concentrated mass elements. The 
spring tines were modeled with plate elements. A 
torsional spring was added at the pivot point to bet- 
ter represent the experimentally determined pitch 
frequency. 


Clamping 

block 

(solid elements) 


Plunge spring tines 
- (plate elements) 

( \ 


Wing 

(solid elements, 
concentrated 
masses) 



Pitch axis 
(rigid elements) 


Elements where 
piezoelectric 
ceramics are 
located 


Pitch spring tines 
(plate elements) 


Figure 15. Finite element model. 


Observation of the physical system in motion in- 
dicated that the boundary condition for the plunge 
spring tines was clamped or guided. The finite el- 
ement model closely predicted the plunge frequency 
when this boundary condition was enforced. Further 
observation of the motion indicated that the pitch 
spring boundary condition at the pitch pivot point 
was stiffer than a cantilever, while the other end of 
the tine looked cantilevered (ref. 21). The exagger- 
ated sketch of figure 16 depicts this phenomenon. 
From experimental frequencies and calculated iner- 
tias, spring stiffness constants were computed for 
three wing configurations: the primary wing struc- 
ture alone, one with the balsa wood wing extension, 


and one also containing a 0.007-lbm mass ballast. 
From the average stiffness value, the cantilevered 
stiffness of the existing finite (dement, model was 
subtracted. The resulting stiffness was included in 
the finite element model by means of an explicitly 
modeled spring at the pitch pivot point. 


Undeformed 



Figure 16. Sketch of boundary condition influence's on 
deformation. 

The locations on the spring tines where piezo- 
electric plates were bonded were defined as compos- 
ite plates, where the layers of piezoelectric ceramic 
were modeled as temperature-dependent materials. 
The parallel constitutive relations of thcrmomechan- 
ical and electromechanical systems allowed the volt- 
age applied to the ceramics to be represented by an 
applied temperature field. 

Aerodynamic model . Unsteady aerodynam- 
ics were calculated with the doublet-lattice method 
(ref. 22) as implemented in the Aeroelastic Vehi- 
cle Analysis (AVA) system of computer codes. The 
doublet-lattice method is a panel method for solving 
the integral equation relating the normal wash and 
the aerodynamic loading for lifting surfaces in sub- 
sonic flow. Discrete lifting elements, consisting of an 
oscillatory doublet line and a horseshoe vortex, ap- 
proximate the loading. The steady- flow effects are 
represented by the vortex; the doublet represents the 
incremental effects of oscillatory, unsteady motion. 
AVA uses the modal displacement vectors to cal- 
culate the generalized aerodynamic forces (GAF’s) 
at discrete reduced frequencies. The program out- 
put is a table for each reduced frequency, where the 
columns of the table correspond to modal and con- 
trol deflections, while the rows correspond to modal 
pressures or forces. Because the wing is rigid over the 
airspeed range of interest, the modal displacements 
are input at six points along the leading and trailing 
edges of the primary wing section. The aerodynamic 
model, shown in figure 17, has 5 chordwisc boxes and 
10 spanwise boxes for a total of 50. The GAF's were 
calculated at Mach 0.05 for eight, values of reduced 
frequency ranging from 0.001 to 2.0. 
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Figure 17. Doublet-lattice aerodynamic box layout. 


Rational function approximations for the 
aerodynamics . The aerodynamics produced by the 
doublet-lattice code are transcendental functions of 
reduced frequency. In order to incorporate them into 
the state space equations of motion, they must be 
approximated by rational functions of the Laplace 
variable s (ref. 23). Equation (B59) illustrates the 
second order approximations made. The Integration 
of Structures, Aerodynamics, and Controls (ISAC) 
system of codes (ref. 24) was used to perform these 
approximations and generate the resultant s - plane 
GAF’s. The objective of the fit is to determine the 
coefficients such that the approximation best fits the 
tabular data in a least-squares sense subject to a 
set of linear equality constraints imposed upon the 
coefficients. The constraints imposed for this model 
are that the approximations will exactly match the 
tabular values at zero reduced frequency for each of 
the modes. Figure 18 shows the GAF’s plotted as a 
function of reduced frequency and the results of the 
approximation. 

Finite difference program for generating 
the actuating matrix. The actuating matrix T is 
calculated by a finite difference program. This code 
calculates the second derivative at the center of each 
of the structural elements by using the displacements 
two node points from both sides of the element. 
The only elements included in this calculation are 
those laminated with the piezoelectric plates. The 
actuating plates near the root of a plunge spring 
consist of six elements, each 0.25 in. long and 1 in. 
wide. The displacements are assumed to be constant 


across the width of the spring tine, so only one 
row of displacements along the length is used in the 
calculations. 

Finite difference program for modeling the 
strain gauge. Strain gauges configured to measure 
cantilever bending arc also governed by the behavior 
of the second derivative of the motion, taken with 
respect to the lengthwise coordinate (ref. 25). The 
same basic program used in calculating the actuator 
influence matrix was modified to calculate the strain 
gauge coefficients. The second derivative was calcu- 
lated at locations near the root of the plunge spring 
and multiplied by the spring tine thickness to predict 
the strain on the surface. 

Generating the state space equations of 
motion. The equations of motion were assembled 
with MATRIXx, a commercially available software 
package from Integrated Systems, Inc. (ref. 26). The 
continuous, open loop model was generated in first 
order form. The procedure, given in table IV, shows 
the details of the A and B matrix calculations and 
the C matrix for strain gauge measurements. These 
matrices arc then discretized with the appropriate 
sample rate. 

Modeling the Control Computer 

Dynamics 

The influences of the zero order hold and one 
time step delay on the closed loop system were ex- 
amined using Matlab, a commercially available soft- 
ware package from Math Works, Inc. (ref. 27). The 
influence of the sample rate was also determined. 
MATRIXx implicitly models both the zero order 
hold and a one-sample delay automatically when the 
discretization command is used. 

Illustrated in figure 19 for a gain feedback control 
law, the control computer introduces its own dynam- 
ics into the feedback path. The digital controller im- 
plementation scheme shifts the output data by one 
sample and applies a zero order hold. The frequency 
response of the digital controller is different for differ- 
ent sample rates. Using a sample rate that emulates 
analog derivative feedback near the frequency of con- 
trol interest allows the system to simulate derivative 
feedback, despite having only displacement measure- 
ments (ref. 28). The current control law utilizes the 
dynamics of the implementation scheme, requiring 
only gain feedback. The frequencies of concern lie 
between 7.9 and 11.1 Hz. Figure 19 shows that for 
a 20-Hz sample rate, the phase is —270°, or +90° 
at 10 Hz. Thus, the phase characteristics simu- 
late a derivative in the frequency range of interest, 
near 10 Hz. 
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Figure 18. Rational function approximations to the generalized aerodynamic forces. 



Table IV. MATRIXx Command File for Assembling 
the Equations of Motion 

// 

//eom.mat 

// 

// Procedure file used to generate open 
loop equations of motion 
// 

// Gamma is the finite difference program 
output matrix which 
// calculates the piezoelectric 
structural influence 
// 

// Cstrain is the second derivative of 
each of the mode shapes 

// with respect to lengthwise coordinate 

// 

// Physical parameters 

chord=4 . 0 

rho=. 11468e-6 

d31=-6 . 35e— 9 

qbar= . 5*rho*v**2 

// 

// Structural stiffness, damping, and mass 

// 

omega=[7.8;10.S9] 
f req=2* pi* omega 
dsi= [ .017 .055] 
ms= [eye (2) ] 
ks=diag(f req. **2 ) *ms 
ds=diag(dsi) *diag(f req) 

// 

// Aerodynamic stiffness, damping, and mass 

// 

ma=qbar* (chord/2/v) **2*a2 
da=qbar* chord/2/v* al 
ka=qbar*a0 
// 

// Combining aerodynamic and structural 
matrices 
// 

m=ms+ma 

minv=inv(m) 

d=ds+da 

k=ks+ka 

// 

// Assembling the state space matrices 

// 

a= [0* one s ( 2 , 2 ) eye (2,2) ; -l*minv*k-l*minv*d] 
// 

bs=— 1 M3 l*minv* gamma 
b=1.2*[0;0;bs] 
zbar= . 008 

c=— 1.3*zbar*cstrain 
s= [a b;c 0*ones(2,l)] 



Frequency, Hz 

Figure 19. Frequency response of implicit controller dynamics 

for 20-Hz sample rate (time delay and zero order hold). 

Scaling and Correction Factors 

The experimental setup contains amplifiers, dis- 
cretizations, etc., which must be included in the an- 
alytical model if a controller design is to be applied 
to the physical system. The strain produced on the 
model is measured by gauges that produce voltages. 
These voltage levels are insufficient for the digital 
controller to discern. Thus, an amplifier with a gain 
of 100 is introduced into the strain path. It must 
also be kept in mind that the strain is not actually 
fed back, but a voltage proportional to strain. Any 
control law generated must account for this factor. 
The feedforward path, from the control computer to 
the piezoelectrics, also contains an amplifier, which 
multiplies the input by 25. This gain is included in 
the computation of control matrix B. 

The equations generated contain errors due to 
inability of theory to predict physical phenomena, 
shortcomings in methods used, neglected terms, and 
nonidealities of the physical model. The computa- 
tion of the strain gauge values by a finite difference 
technique had an error of 30 percent at zero fre- 
quency. Because the differencing was performed on 
plate elements very near the clamped boundary con- 
dition, it was determined that the values yielded were 
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inaccurate. The strain gauge equations were scaled 
by 1.3 to account for this. Zero frequency gains were 
also computed experimentally for the transfer func- 
tion from the piezoelectric voltage to the strain gauge 
output. They were off by 20 percent. This was antic- 
ipated because of the unmodeled bonding layer. The 
control matrix was multiplied by 1.2 in an attempt 
to correct for this difference. 

Analyses 

Analyses were performed with the aeroservoelas- 
tic equations of motion derived and implemented as 
detailed in the previous section. Results from these 
analyses are discussed in three sections: (1) finite el- 
ement model studies, (2) aeroelastic analyses of the 
open loop system, and (3) control law design and 
closed loop analysis. 

Results from system identification tests, which 
will be detailed in the experimental section of this 
paper, were incorporated into the aeroservoclastic 
equations of motion prior to the open loop analysis. 
These experimental results include natural frequen- 
cies, structural dampings, and scale factors for the 
actuator and the sensors. 

Finite Element Analysis 

The finite element model served several purposes: 
(1) the structural matrices were generated, (2) nat- 
ural frequencies and mode shapes were calculated 
by performing a normal modes analysis, (3) a para- 
metric study was performed to design the mass bal- 
last, and (4) parameter variations were performed 
to determine the placement of the actuating plates 
necessary to obtain the maximum control effect. 

After the construction of the structural mass and 
stiffness matrices, normal mode analysis was per- 
formed to generate natural frequencies and mode 
shapes: 

[K, -A,M,] {^} = 0 (4) 

The natural frequencies were calculated from 

f t (Hz)=^sZK (5) 

The resulting undamped mode shapes {\Fj} are or- 
thogonal and, like any eigenvectors, can be arbitrar- 
ily scaled. It is a common practice in aeroelastic 
modeling to scale them such that a unit generalized 
structural mass matrix is generated: 

M. s = I*]” 1 M, [*] = [/] (6) 


Predictions of the natural frequencies are pre- 
sented in table V for various pitch boundary con- 
ditions with and without piezoelectric plates incor- 
porated. The pitch spring boundary condition was 
enforced first as a cantilever and second as a com- 
bination of cantilevered and guided boundary condi- 
tions. The combination of boundary conditions was 
enforced to increase the pitch spring stiffness to the 
value indicated by experimental results. The finite 
element model was augmented with the piezoelectric 
actuators and reanalyzed. The presence of the piezo- 
electric plates increased the open loop stiffness; the 
plunge frequency increased by 15 percent. 

Bernoulli-Euler beam theory calculations pro- 
vided valuable benchmarks for the frequency results. 
Reasonableness checks were performed by compar- 
ing idealized beam theory for cantilevered and guided 
boundaries with the finite element results. The in- 
fluence of the piezoelectric elements on the mass 
and stiffness properties is neglected in these calcula- 
tions. The two degrees of freedom are considered as 
completely decoupled. The plunge spring was mod- 
eled with a cantilevered boundary condition, while 
the pitch spring stiffness was computed for canti- 
levered and then guided boundary conditions. The 
frequency predictions are presented in table V. 

The most accurate system model was deemed to 
be the finite element representation that included the 
combined guided and cantilevered pitch spring model 
with the piezoelectric plates in place. The vibration 
mode shapes are shown in figure 20. The first 
mode, designated plunge because of the dominance 
of translational motion, was predicted at a frequency 
of 7.8 Hz. The second mode, which is characterized 
by the pitching of the wing relative to the mount 
system, has a natural frequency of 10.9 Hz. 




Figure 20. Finite element analysis vibration mode shapes. 
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Iable V . Analytical Predictions of Natural Frequencies 
With and Without Piezoelectric Actuators 



Pitch 

spring 

boundary 

condition 

Plunge 

mode 

frequency, 

Hz 

Pitch 

mode 

frequency, 

Hz 

Without, piezoelectric plates: 




Beam theory 

Cantilevered 

7.2 

6.7 

Beam theory 

Guided 

7.2 

13,1 

Finite element model 

Cantilevered 

6.6 

9.2 

Finite element model 

Combined 

6.8 

10.6 

With piezoelectric plates: 




Finite element model 

Combined 

7.8 

10.9 


Prior to the analysis presented above, the mass 
and inertia properties of the physical configuration 
had to be settled. A study to determine the amount 
of mass ballast desired was performed with the fi- 
nite element model. The natural frequencies of the 
system are dependent upon the mass distribution: 



The total mass of all components supported by the 
plunge spring tines, m, is used in the calculation of 
the plunge frequency. The inertia used in the cal- 
culation of the pitch frequency is the inertia of the 
entire wing about the pitch axis. The mass ballast, 
located at the trailing edge, contributes significantly 
to the inertia. The flutter speed is shown in appen- 
dix A to be sensitive to frequency separation. Thus, 
the flutter characteristics can be adjusted by slight 
modifications to the mass ballast. Table VI compares 
the analytical predictions of the natural frequencies 
for mass ballasts from 0.005 lbm to 0.011 lbm. It 
also shows the natural frequencies of the primary alu- 
minum wing with no extension or ballast and of the 
wing with only the balsa w r ood wfing extension. 

Another study was performed to determine the 
optimal placement of the actuating plates. By us- 
ing 1-in-long piezoelectric segments, a pair of actua- 
tors was analytically placed at different locations on 
the plunge spring tines. Table VII shows the dis- 
placement generated at the u free” end for different 
actuator locations and identical applied voltages. In 


Tabic VI. Influence of Mass Ballast on Analytical 
Predictions of Natural Frequencies 



Pitch 

frequency, 

Hz 

Plunge 

frequency, 

Hz 

Primary wing 

25.78 

8.96 

With wing 

13.29 

8.38 

extension only 



Mass ballast, lbm: 



0.005 

11.39 

7.99 

.007 

10.86 

7.79 

.009 

10.65 

7.61 

.011 

10.41 

7.40 


Table VII. Results of Study To Determine 
Actuator Placement 


Distance from root 
to edge of piezo, 
in. 

Plunge mode 
frequency, 

Hz 

Pitch mode 
frequency, 

Hz 

Model 

deflection, 

in. 

No piezo present 

7.02 

12.6 

0 

0.25 

7.77 

12.7 

-9.21 x 10~ 5 

1.25 

7.26 

12.6 

—5.69 x 10“ D 

2.25 

7.02 

12.5 

-1.54 x 10~ 5 

3.25 

7.02 

12.5 

2.59 x 10“ 5 

4.25 

7.30 

12.6 

6.77 x 10~ 5 

5.25 

7.75 

12.6 

9.47 x 10~ 5 


reference 7 it was found that the actuators should be 
placed in regions of high strain for the mode that is 
of control interest. By virtue of both maximum de- 
flection and maximum strain, this investigation indi- 
cated that the actuating plates should be placed near 
either end of the clamped spring tine. 
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Open Loop Aeroelastic Analysis 

Flutter analysis of the model was conducted by 
analyzing the open loop aeroelastic equations of mo- 
tion at a given density for various velocities (ref. 29). 
As the velocity changes, the relative influence of the 
aerodynamic and structural contributions to the in- 
ertial, damping, and stiffness characteristics of the 
system changes. The velocity root locus plot of fig- 
ure 21 shows the open loop flutter characteristics. 
The plot traces the roots of the system as the air- 
speed is increased. The horizontal axis is the real 
part, while the vertical axis is the imaginary part. 
The imaginary axis represents the point of neutral 
stability or zero damping, where, theoretically, re- 
sponses will neither converge nor diverge. Flutter, 
defined as an oscillatory divergence, is represented 
on a root locus plot by an eigenvalue crossing this 
axis into the right half plane. 



Real part, 1/sec 


Figure 21. Root locus (eigenvalue locations as velocity is 

varied). 

Figure 21 predicts the behavior of the plunge and 
pitch modes for sea level density. The frequencies 
of the two modes migrate toward one another as 
the aerodynamics couple the two modes. When the 
frequencies are close together, the modes interact 
with one another and the system is driven unstable, 
shown by the plunge mode eigenvalue crossing into 
the right half plane. The predicted flutter mechanism 
involves the coalescence of the plunge and pitch 
modes at a velocity of 560 in/sec and at a frequency 
of 9.1 Hz. 


An alternate method of expressing the same data 
is to plot the frequencies and the damping ratio as 
functions of velocity. This method will be discussed 
later in the comparison of analysis and experiment. 

Control Law Design and Closed Loop 

Analysis 

The control law design requirement for the cur- 
rent study was to suppress flutter over a velocity 
range. There was no performance requirement be- 
sides stability. Because of the simplicity of both the 
design objective and the test article, the aeroelastic 
phenomenon was controllable t hrough a single-input 
single- out put control law. Gain feedback (ref. 30) 
utilizing the dynamics of the discretization process 
was employed. The strain-proportional voltage was 
the input to the control law. 

Saturation places limits on the implementable 
feedback gain. In a traditional aircraft control 
scheme, the gains are restricted by aerodynamic sat- 
uration of the control surfaces as they stall or by the 
limits of hydraulic actuators. In this experiment, the 
piezoelectric actuators were capable of handling more 
voltage than the operational amplifier was capable of 
producing. Thus, the maximum voltage of the am- 
plifier, not the limiting voltage of the piezoelectric 
plates, became the limiting value. The amplifier had 
an output limit of 80 V; input voltages were amplified 
by 25. Thus, the maximum input voltage was 3.2 V. 

Maximum gain x Maximum strain response < 3.2 V (9) 

The open loop strain response had a maximum 
measured voltage of 0.097 V for wind tunnel condi- 
tions just below flutter. Thus, the usable gain was 
limited to 33. 

In the implementation of the control law, the sig- 
nal was discretized by a 20-Hz sampler and multiplied 
by the feedback gain. Digital- to-analog converters 
held the output data until the end of the sample pe- 
riod. The output signal was then updated. This 
value was held until the end of the next sample pe- 
riod, when new output data were available. Refer- 
ences 31 and 32 provide more detailed explanations 
of this digital- to- analog conversion, called a zero or- 
der hold. Detailed characterization of the dynamics 
of the control law computer was discussed under the 
analytical modeling section; their influence has been 
included in the following analyses. 

Control law design is traditionally performed in 
the continuous domain. Because the control law com- 
puter dynamics were an integral part of this design. 
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Imaginary part 



Figure 22. Gain root locus of the discrete system at open loop flutter; 580 in/sec with gains from 0 to 120. 


however, the discrete domain model was utilized in 
this investigation. On a diagram of discrete system 
eigenvalues, the stability condition corresponds to 
the location of the roots relative to the unit circle. 
Roots outside the unit circle correspond to instabil- 
ities. That is, the imaginary axis of the continuous 
complex plane maps to the unit circle of the discrete 
complex plane. 

Design models were constructed from the aero- 
servoelastic equations of motion developed in the 
previous section. Because these equations contain 
velocity-dependent terms, models were created at 
several distinct velocities. The equations represent- 
ing the system at the predicted open loop flutter 
condition, determined to be 580 in/sec, were the ini- 
tial design model. The continuous model was dis- 
cretized with a 20-Hz sample rate. A gain root lo- 
cus, constructed by varying the gain from 0 to 120, 
is shown in figure 22. Each eigenvalue trace be- 
gins at the open loop system values that corre- 
spond to a feedback gain of zero. One pair of roots 
shown in figure 22 is unstable for the open loop 
case, as this velocity corresponds to the open loop 
flutter condition. The flutter mode eigenvalues sta- 
bilize for small feedback gains, since they migrate 
inside the unit circle almost immediately. As the 


gain increased, the eigenvalues continued to migrate 
within the unit circle for feedback gains up to 108, 
where one destabilized again. 

The stability criterion can be expressed as a limit 
on the magnitude of the eigenvalues. The magni- 
tude of the largest eigenvalue must be less than 1.0 
for the system to be stable. Figure 23 shows the 
value of the maximum magnitude of the eigenvalues 
plotted against feedback gain. The design model, 
linearized at 580 in/sec, stabilized for gains higher 
than 14 but less than 108. As the gain increased 
from 0 to approximately 45, the system became 
more stable; additional gain did not decrease the 
eigenvalue magnitude. This model, however, repre- 
sented the system at only one airspeed. The same 
figure shows the variation with gain for several air- 
speeds. For 1 in/sec airflow, the model was open loop 
stable. As the gain increased, there was very little 
change in the eigenvalue magnitude until the gain 
reached 102. The magnitude increased and the sys- 
tem was driven unstable for gains of 104 and above. 
The maximum airspeed for which gain feedback sta- 
bilized the system was found to be 1300 in/sec. As 
indicated in the figure, gains between 103 and 108 
stabilized the system at this airspeed. If physi- 
cally attainable and no stability margins are required, 
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103 is the optimal gain. Practical limitations on the 
gain, however, did not allow a gain of this magnitude 
to be implemented. 


0 > 

=} 



Figure 23. Maximum complex absolute value as a function of 

feedback gain for various airspeeds. 

A velocity root locus was constructed with the 
largest allowable feedback gain, 33. Figure 24 shows 
the traces of the eigenvalues as the velocity is in- 
creased from 0 to 700. As with the gain root locus, 
a stable system has eigenvalues all lying within the 
unit circle. The flutter mode initially stabilized and 
then slowly began to migrate back toward the unit 
circle. The root crossed the stability boundary, pre- 
dicting closed loop flutter at 648 in/sec. The pre- 
vious graph (fig. 23) showed the variation of maxi- 
mum eigenvalue magnitude as a function of gain for 
this velocity. The influence of increasing velocity can 
also be seen on this graph by examining the dashed 
vertical line representing a feedback gain of 33. At 
1 in/sec, the eigenvalue is just below the stability 
point, 1.0. The next higher velocity plotted is the 
open loop flutter speed, 580 in/sec. The eigenvalue 
magnitude has decreased, indicating the stabilizing 
effect of the feedback. For a velocity of 648 in/sec, 
the trace intersects the stability boundary. Recall 
that this is the closed loop flutter speed. For veloci- 
ties greater than the closed loop flutter speed, there 
is a substantially increased magnitude of instability. 
This is illustrated by the data corresponding to a 
velocity of 1300 in/sec. 

Initial parametric studies, presented in appen- 
dix A, indicated that separation of the zero-airspeed 



Figure 24. Discrete system root locus for varying velocity of 
the closed loop system with feedback gain of 33. 

natural frequencies, or undamped natural frequen- 
cies, would have the effect of delaying the onset of 
flutter. A comparison between the open and closed 
loop eigenvalue magnitudes for increasing velocities 
is presented in figure 25. The effect of the gain feed- 
back is shown to separate the natural frequencies as 
the aerodynamic influence grows (i.e., the velocity 
gets larger, increasing the magnitude of the aero- 
dynamic contributions to the mass, damping, and 
stiffness matrices). 


Closed loop flutter 



Figure 25. Open and closed loop complex absolute value of 
maximum system eigenvalue versus velocity with 20-Hz 
sample rate and feedback gain of 33. 


Analytical Results Summary 

The finite element results indicated that the 
piezoelectric plates should be placed near the ends 
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of the spring tines. The natural frequencies pre- 
dicted were 7.8 and 10.9 Hz. The results of the open 
loop aeroelastic analysis indicate the onset of flutter 
at 560 in/sec. With gain feedback at a gain of 33, and 
control law computer dynamics imparted by a 20-Hz 
sample rate, the closed loop flutter speed is predicted 
to improve by 15.7 percent. 

Experiments 

Experimental results from bench tests and wind 
tunnel tests are presented in this section. Bench tests 
yielded scaling factors for the feedback sensor and ac- 
tuator and provided zero-airspeed system identifica- 
tion results. The wind tunnel results to be discussed 
include system identification and flutter testing. 

Experiments were performed at various points 
during this project. Static testing was performed on 
the model to check and determine gains within the 
open loop system. This was a very useful test because 
of the many component parts that were necessary to 
conduct the sensing and actuating. Open and closed 
loop flutter tests were conducted and the results 
compared with one another as well as with analytical 
predictions. 

Static Testing 

Static testing was performed on the model to 
determine open loop system gains. To validate and 
correct the mathematical model, an experiment was 
devised to check the strain gauge coefficients. A 
known displacement was applied to the clamping 
block (fig. 9); the strain was measured. Applying the 
same amount of displacement to the mathematical 
model yielded a strain 30 percent smaller. The 
sensor equations associated with the strain gauge 
were increased to give the correct dc value. With 
the strain equation yielding the experimental value, 
a constant voltage was applied to the piezoelectrics 
and a strain was measured. The mathematical model 
predicted 20 percent less strain than the measured 
value. The control matrix was then scaled to yield 
the correct value. 

System Identification Testing 

System identification testing w r as performed us- 
ing several techniques to extract modal frequencies, 
dampings, transfer functions, and general system be- 
havior. There were several experimental and ana- 
lytical techniques used to extract these system pa- 
rameters (fig. 26). Impulse response functions of the 
accelerometer generated by hammer tests proved to 
be the most reliable means of extracting the natu- 
ral frequencies of the system, but could be used only 



Random 
excitation 
voltages 
applied to 
piezoelectric 
plates 



Figure 26. Experimental system identification techniques. 


at zero airspeed because the wing could not be di- 
rectly accessed while in the tunnel. Additionally, the 
amount of disturbance introduced to the flow by the 
presence of the accelerometer and its lead wire dras- 
tically altered the aerodynamic behavior. A second 
technique, employed to obtain a more dramatic re- 
sponse from the pitch mode, was to pluck the spring 
tine and record the free-decay data. This free-decay 
technique was effective at low airspeeds, where the 
plucking did not perturb the model enough to induce 
large oscillations. The third technique was to excite 
the model by applying random voltage to the piezo- 
electric actuators. There were several advantages to 
this method — most importantly, the actuator influ- 
ence was included in the results. Also, because the 
amplitude of the input could be carefully controlled, 
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this technique could be used throughout the test en- 
velope without inducing flutter. By examining the 
recorded input and recorded sensor output, transfer 
functions could be calculated. 

Two of the three methods involve the calculation 
of transfer functions. This is accomplished by read- 
ing the time histories into Matlab and transforming 
them, by fast Fourier transforms (FFT), into the fre- 
quency domain. The autospectrum of the input <£ uu 
and the cross-spectrum of the output with the input 
&uy are calculated: 


<fc uu = FFT (u) x FFT (u) 

(10) 

$ uy = FFT (u) x FFT (y) 

(11) 


The transfer function is the ratio of the cross- 
spectrum over the autospectrum: 


Details of the specific tests are now presented. 

Impulse tests. Hammer taps to produce impulse 
inputs were used at zero airspeed to extract the 
system natural frequencies for various configurations 
of the model prior to its being mounted in the wind 
tunnel. 

The first set of impulse tests was designed to 
extract the uncoupled pitch mode frequency; thus, 
the plunge degree of freedom was constrained. The 
model was configured with and without the balsa 
wood extension and with various amounts of mass 
ballast during these tests. Data were taken for 16 sec 
at 64 samples/sec. Overlap averaging of several runs 
was performed to obtain cleaner data. The results 
of these experiments are given in table VIII and 
compared with the analytical predictions based on 
the finite element model. These results were utilized 
to verify the mass and stiffness predictions obtained 
from the finite element model. 

Table VIII. Influence of Mass Ballast on the Pitch Frequency 



Analysis 

Experiment 

(uncoupled) 

Primary wing 

25.78 


With wing 

13.29 

12.8 

extension only 



Mass ballast, lbm: 



0.005 

11.39 

10.4 

.007 

10.86 

9.4 

.009 

10.65 

No data 

.011 

10.41 

8.25 


A second set of experiments was performed uti- 
lizing the impulsive input, which allowed motion 
in both the plunge and pitch degrees of freedom. 
Figure 27 shows time histories of hammer input, ac- 
celerometer response, and strain gauge output. The 
power spectral density of the acceleration response 
(fig. 28(a)) indicates that the natural frequencies of 
the final configuration are 7.9 and 11.1 Hz. Table IX 
compares the natural frequencies before and after 
the actuator elements were added to the finite ele- 
ment model with the experimentally determined val- 
ues. The structural damping of the plunge mode was 
also determined from these data by taking the ratio 
of the frequency width of the peak at the half ampli- 
tude and the natural frequency. The damping ratio 
is half of this value, 0.017. The power spectrum of 
the strain response (fig. 28(b)) does not define the 
modes as well; however, the natural frequencies can 
still be observed from these data. 

Free- decay tests . Free-decay tests were used to 
extract the open loop damping of the pitch mode and 
to compare the open and closed loop dampings of the 
plunge mode. 

The first set of tests was conducted by constrain- 
ing the model in plunge and then plucking the pitch 
spring tine. This method was utilized because the 
impulse response testing failed to extract data suit- 
able for determining the pitch mode damping. Ran- 
dom input tests could not be used, as there was no 
actuator for this degree of freedom. The usefulness 
of this test was limited to determination of the pitch 
mode damping and frequency because the input sig- 
nal cannot be recorded. Thus no transfer functions 
between the input and output can be derived with 
this method. A resultant time history from the free- 
decay testing is presented in figure 29. Accelera- 
tion response of the open loop system was generated 
by applying seven impulses during 18 sec. Each re- 
sponse was fully decayed before the next was applied. 
These data were analyzed with the logarithmic decre- 
ment technique, which is detailed in appendix E. 
The damping ratio which is half of the structural 
damping g , was determined to be 0.055. 

A similar test was performed on the open and 
closed loop systems at zero airspeed without the 
plunge mode being constrained. Although the pitch 
mode was not constrained either, the model was 
perturbed by plucking the model at the clamping 
block, thus primarily exciting the plunge mode. The 
free-decay results showed damping improvements due 
to the controller. The strain responses to pluck tests 
of the open and closed loop systems were compared in 
figure 30. Both responses were normalized such that 
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(a) Hammer input. 



(b) Accelerometer output. 



(c) Strain gauge output. 

Figure 27. Time histories from hammer tap system identifi- 
cation test. 
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(a) Normalized power spectral density of acceleration due to 
hammer. 



(b) Normalized power spectral density of strain due to 
hammer. 

Figure 28. Frequency domain analysis of hammer test transfer 
functions. 


Table IX. Natural Frequencies With and Without 
Piezoelectric Actuators — Analytical 
and Experimental Results 



Plunge mode 
frequency, 
Hz 

Pitch mode 
frequency, 

Hz 

Without piezoelectric plates: 



Finite element model 

6.8 

10.6 

Ground vibration test 

6.9 

12.3 

With piezoelectric plates: 



Finite element model 

7.8 

10.9 

System identification 

7.9 

11.1 
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Figure 29. Time history of acceleration response to plucks at 
zero airspeed (model constrained in plunge). 



Figure 30. Open and closed loop strain response to free-decay 

test at zero airspeed. 

the magnitude of the first peak was 1.0; the data 
were obtained by using 20-Hz sample rates. From 
this plot, the structural damping is shown to have 
been increased by the presence of the controller. 

Random input tests . The open loop system 
cannot be identified fully with either of the tech- 
niques described above because they do not include 
the influence of the piezoelectric actuator. To obtain 
system transfer functions, the actuators to be used 
for control must also be used in the system identifica- 
tion. This criterion was satisfied by tests that input a 
random signal to the piezoelectric plates. The input 
signal was random white noise with a Gaussian dis- 
tribution and a zero mean value. The amplitude of 
the signal, limited to 3.2 by the operational amplifier, 


was adjusted to be as large as possible at the differ- 
ent velocities tested. The larger the signal was, the 
better the data that were obtained, in terms of co- 
herence. Large excitations near the flutter velocity, 
however, drove the model unstable. 

Random excitations were used to examine open 
loop behavior as the wind tunnel velocity was in- 
creased. Figure 31 shows the magnitude of the strain 
versus frequency at eight subcritical airspeeds. The 
plots are dominated by the plunge mode. The se- 
quence of pictures shows that as airspeed is initially 
increased, the peak magnitude decreases and the 
width remains fairly constant. Thus, the half magni- 
tude point falls lower on the curve, that is, at a wider 
point on the peak. Qualitatively this indicates that 
the damping in this mode increases. At 394 in/sec, 
the plot shows that the damping is decreasing. The 
response gets progressively less damped as the flut- 
ter speed is approached. This corresponds to trends 
in the analytical model as shown by the root locus 

(fig- 21). 

Open Loop Flutter Testing 

The risks associated with flutter testing are mini- 
mized in this experiment by the unique design of the 
test article. Because the mount system is outside the 
tunnel and the model is small, it is possible to stop 
flutter by manually taking hold of the flexible springs 
or the clamping block. 

The flutter tests were conducted by placing the 
model in the wind tunnel at zero airspeed and then 
increasing the velocity. Data were taken for several 
minutes at various airspeeds. During this time inter- 
val the model sat at the tunnel condition. Flutter 
points were defined as the lowest airspeed at which 
the magnitude of the oscillations diverged within the 
data-taking interval. The turbulence within the tun- 
nel was relied upon to be sufficient to perturb the 
model. The increments in velocity were made smaller 
and data taken more frequently as the speed got near 
the predicted flutter value. Flutter was encountered 
experimentally at 580 in/sec; the frequency of the 
oscillation was 9.4 Hz. A time history of the strain 
gauge during a run in which flutter was encountered 
(fig. 32) shows the divergent oscillations that begin 
growing at 4.5 sec and continue to grow until the 
maximum possible amplitude is reached at 9.0 sec. 
At this amplitude, safety stops of the tunnel inhib- 
ited the motion of the model so that it would not 
be destroyed. Frequency domain analysis, performed 
by taking fast Fourier transforms of the strain time 
history data, indicates a flutter frequency of 9.4 Hz 
(fig. 33). 
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Table X. Control Effort Expended To Suppress 
Flutter for Increasing Velocity 


Velocity, 

in/sec 

Control effort, 
V 2 

580 

122.8 

590 

172.1 

630 

196.9 

670 

228.0 

710 

*4888.3 


* Flutter encountered. 


Figure 32. Open loop flutter point; 580 in/sec. 



Figure 33. Frequency domain analysis at open loop flutter 
point; 580 in /sec. 


Closed Loop Flutter Testing 

The majority of the closed loop flutter testing 
was initiated by activating the control law at zero 
airspeed. Proceeding in the same manner as the 
open loop flutter testing, the speed was increased 
until flutter was encountered. The closed loop flutter 
velocity was 697 in/sec at 9.7 Hz. 

Increasing airspeed required more control energy 
to be exerted. As a quantitative indicator of this, the 
control effort was defined as the sum of the discrete 
control command signals squared. The control effort 
expended at several airspeeds is presented in table X. 
The control effort for this study has been defined 
as the sum of the absolute values squared of the 
control computer output. As the velocity increases, 
the system responses and the control effort increase. 
The data analyzed at 710 in/sec are actually above 
the closed loop flutter speed. This is the cause of 
the dramatic increase in control effort observed. At 
this velocity, there is not enough control authority to 
suppress flutter. 

The closed loop data were compared with open 
loop data. Figure 34 shows time histories of the 


strain response for the two cases. These data were 
obtained just below the open loop flutter speed, at 
approximately 575 in/sec. The magnitude of the 
peaks for the closed loop case is decreased below 
the magnitude of the open loop peaks. Because 
of limitations in the controller programming, the 
controller update rate and data sampling must be 
consistent. The closed loop data were therefore 
obtained at a 20-Hz sample rate. The open loop data 
appear smoother due to a higher sampling rate. 

The results of the flutter experiments and analy- 
ses are summarized in figure 35. The top graph of fig- 
ure 35 shows the natural frequency variation as veloc- 
ity is increased. Examining the zero-airspeed data, 
the measured open loop values (7.9 and 11.1 Hz) 
are slightly higher than the analytical predictions of 
these frequencies (7.8 and 10.9 Hz). The higher fre- 
quency mode is the pitch mode; the lower frequency 
mode is the plunge mode. The analytical values for 
the two frequencies approach each other as the veloc- 
ity increases. When frequencies coalesce, flutter gen- 
erally onsets soon thereafter. The experimental open 
loop flutter frequency, 9.4 Hz, is slightly larger than 
the frequency that the analysis predicted for flut- 
ter. The final point on this graph is the experimental 
closed loop flutter frequency, 9.7 Hz. 

The lower graph of figure 35 shows the damping 
ratio variation with velocity. Flutter is indicated on 
this graph by the damping ratio of one of the modes 
going to zero. The measured values of damping ratio 
(0.017 and 0.055) from the zero- airspeed data were 
utilized in the analytical model. The analysis results 
are plotted over the velocity range until flutter occurs 
(i.e., the lower mode crosses the zero-damping axis) 
at 560 in/sec. This value is 3.5 percent conservative, 
as the experimental flutter speed was 580 in/sec. The 
closed loop data at zero airspeed indicate that the 
control law provided an increase in the damping. 
Only the data for the plunge mode are presented. 
Analytical predictions of the damping behavior are 
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Figure 34. Experimental open and closed loop strain re- 
sponse to wind tunnel turbulence just below flutter ve- 
locity (575 in/sec). 
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Figure 35. Analytical and experimental flutter results. 

presented over the entire velocity range, until closed 
loop flutter is predicted to occur at 648 in/sec. The 


experimental value was 697 in/sec. The predicted 
value is 7.6 percent conservative. With single-input 
single-output gain feedback, the flutter speed was 
predicted to improve by 15.7 percent. The controller 
actually achieved an improvement of 20 percent in 
the flutter velocity. 

Conclusions and Recommendations 

This research effort has resulted in the first ex- 
perimental demonstration of flutter suppression em- 
ploying piezoelectric actuators. A wind tunnel model 
was conceived, designed, fabricated, installed, and 
tested. Structural and aerodynamic models were cre- 
ated; the aeroservoelastic equations of motion were 
derived and analyses performed. A digital control 
law was designed based on a discretized model and 
was implemented. Open and closed loop flutter tests 
were conducted, with excellent correlation achieved 
by analytical predictions. 

A two- degree-of- freedom wind tunnel model con- 
sisting of a rigid wing attached to a flexible mount 
system was designed based on preliminary flutter 
analyses. The rigid wing, with a primary section 
made of aluminum, was connected to a cantilevered 
spring tine to control the pitching degree of freedom, 
and the entire assembly was then connected to a set 
of spring tines to control the plunging motion. The 
configuration and dimensions of the model were de- 
signed such that it would flutter well within the op- 
erating envelope of the tunnel, could be safely tested 
within the available test section, and would have sur- 
faces suitable for mounting the piezoelectric plates in 
a bimorph configuration. 

Analytical modeling of the wind tunnel model 
resulted in aeroservoelastic equations of motion. 
The equations were derived from Lagrange’s energy 
method and utilized modal analysis of a discretized 
structural model. The natural frequencies were pre- 
dicted to be 7.8 Hz for the plunge mode and 10.9 Hz 
for the pitch mode. Generalized aerodynamic forces 
were generated by the doublet-lattice method and 
approximated with rational functions. Expressions 
for the generalized forces associated with the control 
inputs were derived based on classic laminated plate 
theory and calculated with finite differencing tech- 
niques applied to the discretized structural model 
results. 

Aeroelastic analysis of the open loop system gave 
a flutter prediction of 560 in/sec. Utilizing the 
implicit dynamics of the control law computer, a gain 
feedback control system was designed with strain as 
the feedback signal. After the gain was optimized for 
the largest stable velocity range and saturation of the 
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electronic hardware was accounted for, a gain of 33 
was employed. The flutter speed for the closed loop 
system was predicted to he 648 in/sec, a 15.7-percent 
increase. 

Experimental results from several system identifi- 
cation tests determined the natural frequency of the 
plunge mode to be 7.9 Hz and that of the pitch mode 
to be 11.1 Hz. The structural dampings associated 
with these modes were also determined. The open 
loop flutter speed was measured at 580 in/sec. The 
analytical prediction was conservative by 3.5 percent. 
Closed loop flutter testing was performed and a flut- 
ter speed of 697 in/sec was obtained. This represents 
a 20-percent improvement from the open loop case. 
The analytical prediction of closed loop flutter speed 
was conservative by 7.6 percent. 

It is recommended that further research be per- 
formed in the area of controlling the aeroclastic re- 
sponses of a vehicle by utilizing piezoelectric actua- 
tors. A more realistic and complex model needs to 
be designed that incorporates strain-actuating ele- 
ments within the wing design. The concept has been 
proven to work, however, it has not yet been shown 
to be viable in terms of real aircraft. Experiments on 
a larger scale are now called for. 

Flutter suppression is not the only aeroelas- 
tic application that may call for secondary actua- 
tors made of adaptive materials. Load alleviation 


currently performed by aerodynamic control surfaces 
may prove to be an ideal application for localized 
strain actuation. By actuating adaptive material ele- 
ments, local strains could be produced to counter the 
loads induced within the structure during maneuver- 
ing. This has the potential to extend the service life 
of aircraft that traditionally undergo high ^-loading 
and to expand operational limits. 

The concept of an adaptive material mission- 
adaptive wing is worthy of investigation. It would 
eliminate the hydraulic problems encountered on pre- 
vious attempts to create a wing that can be shape 
optimized for various flight conditions. 

The applications for which adaptive material will 
be suitable in the future depend heavily on the re- 
search in the materials area. Ceramics, which were 
used in this investigation, are very fragile. Polymers 
are currently not capable of generating the force lev- 
els required for actuating realistic structures. For 
piezoelectrics to be utilized in aeroclastic applica- 
tions in an arena beyond research, a more resilient 
substance than ceramics, or a means to protect the 
ceramics, needs to be developed. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
October 21, 1092 
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Appendix A 

Parametric Studies for Model Design 

During the early stages of model design, sensi- 
tivity studies were performed to determine the in- 
fluences of several key parameters on the flutter 
characteristics. The baseline structural model was 
formulated in terms of natural frequencies and a mass 
matrix. The mode shapes were assumed to be uncou- 
pled rigid body plunge and pitch modes. The mass 
properties of the wing were specified as those of an 
isotropic steel plate with the pitch pivot point at the 
midchord. The mount system stiffnesses for plunging 
and pitching degrees of freedom were initially chosen 
such that the natural frequencies were at 9 and 18 Hz, 
respectively. The baseline aerodynamic model was a 
rectangular planform with a chord length of 2 in. and 
a span of 4 in. 

The parametric studies were performed on the 
baseline model by changing the design variables inde- 
pendently and noting the changes in flutter velocity 
and frequency. Figures Al through A 4 show the re- 
sulting trends for variations of the static unbalance, 
frequency ratio, pitch pivot location, and structural 
damping. The values on the plots for variation with 
static unbalance, pitch pivot location, and structural 
damping have been normalized independently of the 
others. The flutter velocities on a given plot are 
divided by the flutter velocity corresponding to the 
lowest value of the parameter being varied. The flut- 
ter frequencies were normalized in the same manner. 
Thus, on each plot, the value of the left-most point 
will be 1.0 and apply to both the flutter velocity and 
the flutter frequency. 

Static unbalance is defined by the off-diagonal 
terms of the symmetric mass matrix. As the static 
unbalance was increased, the flutter velocity de- 
creased, as shown in figure Al. This is indicative of 
the additional mass coupling increasing the tendency 
for modal coalescence. The frequency trend for the 
same variation showed an increasing frequency. This 
indicates that the higher frequency mode played a 
more significant role in the flutter mechanism as the 
static unbalance term of the mass matrix grew. 

The flutter velocity trend with respect to the ra- 
tio of frequencies is shown in figure A2. The ratio 
of frequencies is defined by the plunge frequency di- 
vided by the pitch frequency. Note that the increas- 
ing ratio represents the distance between the natural 
frequencies decreasing. Two sets of data are plotted 
in figure A2. The first set (1) is for various plung- 
ing frequencies divided by the baseline value of the 
pitching frequency, 18 Hz. The graph shows that as 
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Figure A 1 . Influence of static unbalance on the flutter velocity 
and frequency. 
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Figure A2. Influence of frequency ratio on the flutter velocity. 


— Plunge frequency increasing 

— Pitch frequency decreasing 



the ratio of the frequencies gets larger, the flutter 
speed decreases. The second set (2) of data is for 
the baseline value of the plunging frequency, 9 Hz, 
divided by various values of pitching frequency. The 
trend is also for the flutter velocity to decrease as this 
ratio increases. The relative slopes of the two lines 
indicate that over the range considered, the flutter 
velocity is more sensitive to reductions in the pitch 
frequency than to increases in the plunge frequency. 


The plot of the pitch pivot location (fig. A3) 
indicates that as the pivot point is moved toward 
the trailing edge of the wing, the flutter velocity 
is lowered. Recalling that for these variations the 
center of gravity was located at the midchord, or 1.0, 
locations aft of 1.0 are for statically unstable wings. 
Long before flutter, divergence will have occurred. 
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Figure A3. Influence of pitch pivot location on the flutter 
velocity and frequency. 


Figure A4 shows the change in flutter velocity 
and frequency as the structural damping is increased 
simultaneously in both the plunge and the pitch 
modes. The trends show that the velocity increases, 
as expected, since the eigenvalues at zero airspeed 
will be farther from the instability point (i.e., the 
modes are more stable). The trend is not linear, but 
for low values, each percent of additional structural 
damping raises the velocity by roughly a percent. The 


frequency of flutter is shown to decrease as structural 
damping is added. The addition of structural damp- 
ing lowers values of damped frequencies and brings 
them closer together, predictably lowering the flutter 
frequency. 
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Figure A4. Influence of structural damping on the flutter 
velocity and frequency. 


Changes to the design based on the trends seen 
in these parametric variations greatly reduced the 
number of iterations required in arriving at a final 
configuration. 
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Appendix B 


Equations of Motion for Aeroelastic Systems Utilizing Laminated Plate Theory 

This appendix describes the derivation of Lagrange’s equations of motion for an aeroelastic system controlled 
by piezoelectric actuators. The equations derived in appendix D are utilized for the calculation of the actuator 
components of the equations. The constitutive relations developed in appendix C are also incorporated into 
the potential energy equation. 

Lagrange’s equations of motion require the derivation of expressions for potential and kinetic energy as 
functions of generalized coordinates. The aeroelastic modeling was performed with the orthogonal or undamped 
modes of vibration as the generalized coordinates. 

Plate displacements may be expressed in terms of these coordinates: 

L 

u (x, y, *,*) = £ xi (x it Vi, Zi) q xi ( t ) (Bl) 

i — 1 

M 

V (•£. Vi t) — ^ ^ yi {%ii Vii %i) Qyi (t) (B2) 

i=l 

N 

m (%, V, Z, t) = ^ ^ ^ zi {%ii Vii Z{) Qzi (f) (B3) 

z=l 

where &( Xi y tZ )i are the mode shape vectors in the x-, y-. and ^-directions, respectively, and qt x>yig \i(t) are time- 
dependent generalized coordinates. The ^-translation, u, of the system at the (xo, 2/0: z o) location is expressed 
as 

u (^o, 2/0. zq) = ^xi (^o. yo , zo) Qx i + * x2 (^o. yo. z 0 ) <7x2 + ■ ■ ■ (B 4 ) 

In matrix notation, the translational degrees of freedom are written 
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The potential energy of the system is defined as 
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The constitutive equations, detailed in appendix C, provide the relationship between the stress and the strain: 

a = G (e - A) (B9) 
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The total strain e can be calculated based on midplane strain €q and curvature rc. 

The longitudinal and lateral 

strains and the shear strain are defined as 

du 

£xx — 

ox 

(BIO) 

dv 

£ yy = dy 

(Bll) 

du dv 

lxy =dy + di 

(B12) 

For a plate, the displacements are related to quantities at the midplane (midplane denoted by the subscript 0): 

dwQ 

U = UQ z 

dx 

(B13) 

dwQ 

v = vo z 

dy 

(B14) 

o 

3 

II 

3 

(B15) 

which allows the plate strains to be expressed as 


duo d^WQ 

£xx — ^ Z a 9 

ox ox A 

(B16) 

dvQ d 2 w$ 

yy dy d 2 y 

(B17) 

_ du 0 t dv 0 o _ d 2 w 0 
^ xy dy dx Z dx dy 

(B18) 

or 


£xx = £xx 0 T zk x 

(B19) 

£ yy = E yyo + ZK y 

(B20) 

'Ixy = Ixyo T ZK X y 

(B21) 

where the following definitions are made: 

du o 

£xxo ~ dx 

(B22) 

dv o 

£ yyo- dy 

(B23) 

du 0 dv o 

lxy dy + dx 

(B24) 

d 2 WQ 

KX ~ dx 2 

(B25) 

II 

1 

% 

(B26) 

_ d 2 wo 
Kxy ~ * dx dy 

(B27) 

Strain, expressed as 


£ = £Q + ZK 

(B28) 
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is substituted into the constitutive equations to provide 


<j = G (sq + zk — A) 


The potential energy now becomes 


U 


= ^ JJJ (eo + zk) t G (e 0 + zk — A) dV 


Volume 


U — ~ JJJ (so + zk) T G (eo + zk) dV — ~ JJJ (erg -f zk) t G {A} dV 

Volume Volume 

Note that the midplane strain and the curvature are independent of z. Thus, 

JJJ dv =JI I dzdA 

Volume Area Thickness 

The potential energy can then be calculated: 


°-Uh T 

J G dz 

£0 + £0 T 

J Gzdz 

K + K T 

[ Gzdz 

J 

Area 

- Thickness 


- Thickness 


- Thickness 


£0 


+ k t 

f Gz 2 dz 

J 

KdA ~\ [ [ £ 0 T 

[ GA dz 

+ k T 

[ GzA dz 


- Thickness 

Z J J 
Area 

j 

- Thickness 


J 

- Thickness 


dA 


The integrals through the thickness are defined as: 

-J 


A= I G dz 

Thickness 


/ 


B= I Gzdz 

Thickness 

D = J G z 2 dz 

Thickness 

N a = f GA dz 

Thickness 

= J GzA dz 

Thickness 


The A, B , and D matrices are called extensional, coupling, and bending stiffnesses, respectively. 


(B29) 

(B30) 

(B31) 

(B32) 

(B33) 

(B34) 

(B35) 

(B36) 

(B37) 

(B38) 
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The potential energy, expressed in terms of these integrals, is written 


U = - JJ e 0 T Ae 0 + 6q - h T- k DtzdA ^ sq T k, dT 


Area 


Area 


or 


u 


\SS u r « t >[2 «*>{$}“ 

Area Area 


(B39) 


(B40) 


Referring back to the strain and curvature functions in terms of displacements, an operator V is defined such 
that 

' UQ 


K } “ PI 


V{) 

U?o 


(B41) 


where 


[w = 


d 

m 


n d 

0 ^ 

d d 

By Bx 


0 

0 


0 

0 


0 

0 


0 0 0 - 


B 2 d 2 


fa? By 1 ^ Ox dy J 


d 2 


(B42) 


The general relation from midplane displacements to displacements at a distance z from the midplane is 
easily derived. However, if the displacements are given for the midplane, as is the case for plate element results 
of a finite element model, there is a one-to-one correspondence. The subscripts on the displacement vector are 
thus dropped. 

Recalling the expression for the displacements in terms of generalized coordinates allows the midplane strain 
and curvature to be expressed: 

(B43) 


£ k °} = [P] [*]q(t) 
Substituting this expression into the potential energy expression, 

h / k * 


.T V T 


Area 


A B 
B D 


Ih T * v {Z] 

Area 


dA 


(B44) 


or, moving the generalized coordinate outside the integrals 


where 


- II ^ 


Area 


II 


Area 


1 'T 1 

2 q F2 

(B45) 

B 

[D¥] dA 

(B46) 

D 


N * } 

(B47) 

J 



¥\ is the generalized structural stiffness matrix K s . 
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(B48) 


The kinetic energy development follows a similar path. By definition, kinetic energy T is 

T = \ JJJ p{u} T {u}dV 

Volume 


where 



(B49) 


From equation (B5), the displacement vector is expressed by the mode shapes multiplied by the generalized 
coordinates. The mode shapes arc independent of time, so that the expression for the time derivative of 
displacement is 

{«} = [#] {q} (B50) 

where ( ' ) represents the derivative with respect to time. The kinetic energy is then written 


T = \ JJJ pq 7 ’* T M>qdF 

Volume 


(B51) 


Moving the generalized coordinates outside the integral, the expression becomes 


T = 


^ /// 


pty 1 Vl> dV q 


Volume 


(B52) 


or 

1 r 

T = 2 <T F;jq (B 53 ) 

where 

f 3= JJJ p'& T 'f'dV (B54) 

Volume 

The triple integral over the volume can be expressed as the double integral over the area of the integral over 
the thickness. Assuming that the mode shapes are constant through the thickness, the mass per unit area, 
m 0i is defined as the integral of the density through the thickness. The expression for F;$ is recognized as the 
generalized structural mass matrix: 

M,=F 3 = JJm Q ^ T ^d,A (B55) 

Area 

In formulating the Lagrangian equations of motion, the conservative forces of the system are contained in the 
potential and kinetic energy expressions (ref. 33). The nonconservative forces, namely the aerodynamics, are 
represented by generalized forces Q , : 


d_ fdT\ dT dU dU _ 
dt \d^)~ % + d^ + W t =Ql 

From the previously derived energy expressions, the Lagrange equations become 

„ <9D 

F3q + ~qT- + Fiq=Q; + F2 

which can be transformed to the Laplace domain and rewritten as 

(M.,.s 2 + D s s + K.,) {q/} + q [Q/] {q/} = F 2 


(B56) 


(B57) 


(B58) 
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This is the classical second order aeroelastic equation formulation with the exception of the term on the right- 
hand side; F 2 is the force generated by the actuating strain, which arose from the addition to the potential 
energy expression. It appeared because a mechanical stress was being produced by a nonmechanical strain. 
In the aeroelastic problem, Lagrange’s equations are written as a balance of mechanical energy. Further 
development of the strain actuation matrix is presented in appendix D. 

In order to transform the aerodynamics into the Laplace domain, it was necessary to apply a second order 
rational function approximation to the aerodynamics (ref. 34): 

Q/ - Q/ = [ A() ] / + [ Al l / * (^) + M / s2 (£) (B59) 


Equation (B58) is then expressed 

(Ms 2 + D.s + k) {q/} =] 

where 

K = K„ + g [Ao] 


6=d * + *(£)[H 

m-m s + 9(^) 2 [a 2 ] 

Defining the vectors 

*/, = {q/} 


x h = s M 

these equations can be readily converted to first order form (ref. 35): 


x h 

x h 


0 

-M ! K 


I 

-M ■£> 



(B60) 

(B6I) 

(B62) 

(B63) 

(B64) 

(B65) 

(B66) 


From the development of the strain 
are expressed as: 



actuation matrix in appendix D, the acroscrvoelastic equations of motion 


0 

-M^K 


I 

-M ! D 



0 

-M" l rd :n 



(B67) 


In traditional controls notation (ref. 30), 


x = Ax + Bu 


(BG8) 
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Appendix C 

Development of the Constitutive 
Relations for a Piezoelectric Material 

This appendix develops the constitutive equations 
for a piezoelectric material, taking them from the 
generalized Hooke’s law to a subset utilized in lami- 
nated plate theory and expressed in a form useful in 
the derivation of the potential energy equation. 

Mechanical stresses and strains are related 
through a 6 by 6 compliance matrix in the gener- 
alized Hooke’s law. The constitutive relations for a 
linear elastic material with three mutually perpen- 
dicular planes of elastic symmetry have nine indepen- 
dent entries in the compliance matrix. Classic lami- 
nated plate theory (ref. 35) is used in developing the 
equations, so only the in-plane stresses and strains 
are considered, reducing the independent elements 
to four, as shown in the following equation: 


The nonzero terms for a PZT plate with poling in 
the 3-direction are 


[d] = 


0 

0 


ld 3l 


0 0 0 di5 

0 0 d 24 0 

d 3 2 d 33 0 0 


o- 

0 

0. 


(C8) 


Once again, by using classic laminated plate theory, 
this matrix reduces to a 3 by 3: 


d = 


0 

0 


L^31 


0 d\ 5" 

0 0 
d$2 0 - 


(C9) 


When the electrodes are on the faces of the plate per- 
pendicular to the poling direction, the only voltage 
that can be applied is in the 3-direction. Therefore, 




1 


0 ' 


f 1 


n 

r Exx j 
< £ vv 1 

> — 

-3? 

1 

0 

i 

1 1 
1 ° y 1 

> (Cl) 

E = l 0 > 




1 

G -. 12 - 




1 e 3 . 

Ixy ' 


0 

0 


1 T X y ) 




(CIO) 


which can be expressed as 


Defining A as the actuation strain vector leads to 


£ — S <X (C2) 

The strain vector £ and the stress vector cr are related 
by the compliance matrix S. An alternate method 
of expressing this relationship is through a stiffness 
matrix G: 

cr — Gs (C3) 



f A1 1 

( ^31 ' 

A = < 

A2 

> = d T E = < d 3 2 > 


1a 3 J 

[ 0 . 


(Cll) 


Solving equation (D6) for the stress vector, and 
defining G as the inverse of the compliance matrix, 


The electrical quantities (flux density R and voltage 
per thickness E) are related by the dielectric equa- 
tions through permittivity j3 and impermittivity e 
matrices (ref. 7): where 


where 


E = /3R 

(C4) 






El 

v 2 \E 2 

0 

R = eE 

(C5) 


1 - ^12^21 

1 - V12V21 



G = 

^12^2 

e 2 

0 

V 

(C6) 

1 - v 121^21 

1 - ^ 12^21 

E = — 
t 


0 

0 

G n 


(C12) 


(C13) 


For piezoelectric systems, the coupling of the two 
fields is accomplished by the introduction of the 
strain-charge matrix: 

e = S <r + d T E (C7) 


Note that use of the laminar plate model neglects 
any influence of the d 33 coefficient and that applying 
the voltage parallel to the polarization negates the 
influence of the ^15 coefficient. 
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Appendix D 

Analytical Modeling of the Piezoelectric 
Actuator 

The equations derived in appendix B consist of 
components that are calculated by standard meth- 
ods, with the exception of the force due to the piezo- 
electric actuators F2. This appendix presents a more 
detailed look at the calculation of this matrix, as 
well as the approximations used to implement the 
calculations. 


derivative taken with respect to x will be considered 
negligible, leaving 

[Z>*] r =[o 0 0 0 0] (D 6 ) 

Thus, 

F2= / (■' w) Mh ‘ dA (D7) 

Area 


From equation (B47), the force is given as 


F 2 = 


II [v * f 


Area 



(Dl) 


The rows of the force and moment due to strain 
actuation may be specified 


N A = { 


M\=< 



M A] n 

Ma 2 

Ma 3 


(D 2 ) 


(D3) 


The definitions of the operator and modal matrices 
are recalled from equations (B6) and (B42): 


[ 2 >¥] = 


r £ 0 

0 

d d 

cJy JJx 


0 

0 

0 

_ d 2 
dx 2 

_ d 2 
dy 2 

2 & 

i Ox dy 


[*d 

0 

0 


0 

[*K 

0 


0 1 
0 

[*d 


( 04 ) 

Consider only the out-of-plane displacements (i.e., 
%? Xi and ty yi are zero); the matrix dimension shrinks 
to 


[DH>] T =[000 


d 2 * z 
dx 2 


d 2 * z o a 2 * z 

dy 2 dx dy 


(D5) 


In the design of the flexible mount system, special 
consideration was given to preventing the plunge 
spring tines from deforming in torsion. The piezo- 
electric plates are oriented along the spring tines, 
which deform primarily in the y-direction (i.e., cre- 
ating moments about the x-axis). Therefore any 


The piezoelectrically induced moment was derived in 
equation (B38): 


M a = J GzAdz 

Thickness 

The second row, corresponding to M\ , is 


(D 8 ) 


[t 


^12^12 

~ ^ 12^21 1 - v \2 v 2i 


Thickness 


oj i A 2 1 zdz (D 9 ) 


Recall from the discussion of constitutive relations 
that 


f A1 l 

H 

A 4 = < 

1 

l 0 J 

l 0 ; 


} (hi 


V 3 


(DIO) 


Then 


Ma, 


= [ (.— 2 + 1 ) E 2 d-M 1 

J \ 1 - 1^12^21 / t 


y 3 zdz 


(Dll) 


Thickness 


Defining 


provides 


7= ( 1 1/12 + 1 )E 2 d 3l \ (D12) 

\1 - V\2l^2 1 / 1 


Ma, 


-I 


jV^zdz 


(D13) 


Thickness 


The scalar 7 consists of geometric and material 
properties. For any isotropic lamina, 7 is indepen- 
dent of the location 2, so it may be taken outside 
the integral. The structure is considered to consist 
of several layers composed of different materials. The 
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Bond layer 
0.007. in. thick 


Piezoelectric plate 
0.0075 in. thick 

- — — — — a — 0.015 in. 



Plunge spring tine 
0.016 in. thick 


Figure Dl. Actuator attachment. 


integral over the thickness is decomposed into several 
parts: 


m a . 


(-/ 

\ Layer 1 


Z dz + 72 / z dz 

Layer 2 

+ 73 / z dz + 74 J zdz 

Layer 3 Layer 4 


+ 75 


/ 


\ 


z dz 
Layer 5 / 


v 3 


(D14) 


where the layer geometry is defined in figure Dl. 

The constant 7 is a function of the electro- 
mechanical coupling coefficient d%\. Neither the steel 
nor the bonding compound exhibits any coupling be- 
havior, so for layers 2, 3, and 4, the coupling coef- 
ficients are all zero. As long as oppositely oriented 
voltages are applied to the top and bottom plates, 
the plates are geometrically and electrically similar 
such that the integrals through the layers are equal. 
Therefore, 


and 


F2 = If- ^71 V 3 (2 zt + t 2 ) dA (D17) 

Area 

To implement the actuator equation, numerical inte- 
gration must be performed. The integral equation is 
approximated as a summation over the node points 
of a discretized structural model (ref. 32): 


Area 


„ r n nodes 

\)dA= £ [ 
3 = 1 


(D18) 


where Aj is the surface area associated with each 
node. Therefore 


n nodes / 

F 2 = E 7l(2*t + t 2 ) - 

3 = 1 ' 


dy 2 


AjV 3 (D19) 


There are two modes used in the analysis, so 

/ n nodes 


F 2 = 71 (2zt + t 2 ) < 


nodes / \ 

S (-»/’ 

I n nodes / _ 2 \ 

l ,5 


} Vs (D20) 


/*£+ 1 

Ma 2 = 2 7iV 3 / zdz 

(D15) 

F 2 - I r T ^3iV 3 (D21) 

J Z 


l r 2 J 

After solving the integral, 


where T* is associated with the ith mode: 



Ma 2 = 71 V 3 (2zt + t 2 ) 

(D16) 

n nodes / 9 \ 

r ' -<-•>(. e (-'#)* <D22) 

j — 1 x / 3 
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Appendix E 

Logarithmic Decrement Method for 
Determining Damping 

The logarithmic decrement is defined for a decay- 
ing cyclic system by the ratio of peak magnitudes for 
two cycles, which are n cycles apart: 



(El) 


The relationship of the logarithmic decrement to the 
damping ratio is 



which for small values of damping can be approxi- 
mated 

S = 2ttC (E3) 


giving the formula for damping ratio 




(E4) 


Note that the structural damping g is twice the 
damping ratio C- 
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and tested. The model consisted of a rigid wing and a flexible mount system that permitted a translational 
and a rotational degree of freedom. The model was designed such that flutter was encountered within the 
testing envelope of the wind tunnel. Actuators made of piezoelectric material were affixed to leaf springs of the 
mount system. Command signals, applied to the piezoelectric actuators, exerted control over the damping and 
stiffness properties. A mathematical aeroservoelastic model was constructed by using finite element methods, 
laminated plate theory, and aeroclastic analysis tools. Plant characteristics were determined from this model 
and verified by open loop experimental tests. A flutter suppression control law was designed and implemented 
on a digital control computer. Closed loop flutter testing was conducted. The experimental results represent 
the first time that adaptive materials have been used to actively suppress flutter. They demonstrate that small, 
carefully placed actuating plates can be used effectively to control aeroelastic response. 
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